Shifting Norms - Moving Climate Norms to the New Climatology Period

# load average global climatologies as annual averages
yr_avgs <- str_c(box_paths$oisst, "daily_climatologies/yearly_means/")
clim_82_avg    <- stack(str_c(yr_avgs, "avg_clim_1982to2011.nc "))
clim_91_avg    <- stack(str_c(yr_avgs, "avg_clim_1991to2020.nc"))
clim_shift_avg <- stack(str_c(yr_avgs, "avg_clim_shift_82to91baselines.nc"))

# # load daily_climatologies
# clim_dailys <- str_c(box_paths$oisst, "daily_climatologies/")
# clim_82_daily    <- stack(str_c(clim_dailys, "daily_clims_1982to2011.nc "))
# clim_91_daily    <- stack(str_c(clim_dailys, "daily_clims_1991to2020.nc"))
# clim_shift_daily <- stack(str_c(clim_dailys, "clim_shift_82to91baselines.nc"))


# climatology_lists
grid_list <- list(
  "1982-2011"      = clim_82_avg,
  "1991-2020"      = clim_91_avg,
  "Climate Shift"  = clim_shift_avg)

# rotate them
grid_list <- map(grid_list, raster::rotate)
# convert annual_avgs to stars objects
clim_82_st    <- st_as_stars(grid_list$`1982-2011`)
clim_91_st    <- st_as_stars(grid_list$`1991-2020`) 
clim_shift_st <- st_as_stars(grid_list$`Climate Shift`) 


# Warp to world projection
wgs_84          <- st_crs(4326)
robinson_crs    <- st_crs(54030)
world_moll      <- st_crs("+proj=moll")
clim_82_moll    <- warp_grid_projections(clim_82_st,    projection_crs = "world moll")
clim_91_moll    <- warp_grid_projections(clim_91_st,    projection_crs = "world moll")
clim_shift_moll <- warp_grid_projections(clim_shift_st, projection_crs = "world moll")


# Get temp limits that fit both
overall_max <- max(c( max(clim_82_st[[1]], na.rm = T) , max(clim_91_st[[1]], na.rm = T) ))
overall_min <- max(c( min(clim_82_st[[1]], na.rm = T) , min(clim_91_st[[1]], na.rm = T) ))
temp_limits <- c(overall_min, overall_max)


# Map the two periods
plot_period <- function(in_stars_obj, period_label, crs_projection){
  period_plot <- ggplot() +
    geom_stars(data = in_stars_obj) +
    geom_sf(data = world, fill = "gray90") +
    geom_sf(data = st_transform(world, crs_projection), fill = "gray90") +
    scale_fill_distiller(palette = "RdYlBu", na.value = "transparent", limit = temp_limits) +
    guides("fill" = guide_colorbar(title = expression("Average Annual Sea Surface Temperature"~~degree~C),
                                   title.position = "top", 
                                   title.hjust = 0.5,
                                   barwidth = unit(4, "in"), 
                                   frame.colour = "black", 
                                   ticks.colour = "black")) +
    coord_sf(crs = crs_projection) +
    map_theme +
    labs(title = period_label)
  
  return(period_plot)
}
  

# Plotting the two periods
plot_82 <- plot_period(clim_82_moll, period_label = "1982-2011 Climatology", crs_projection = world_moll)
plot_91 <- plot_period(clim_91_moll, period_label = "1991-2020 Climatology", crs_projection = world_moll)



# Set diverging limits for scale
shift_limits <- max(abs(clim_shift_st[[1]]), na.rm = T) * c(-1,1)

# Plot the shift in the two
plot_shift <- ggplot() +
    geom_stars(data = clim_shift_moll) +
    geom_sf(data = st_transform(world, world_moll), fill = "gray95") +
    scale_fill_distiller(palette = "RdYlBu", na.value = "transparent", limit = shift_limits) +
    guides("fill" = guide_colorbar(title = expression("Sea Surface Temperature Difference"~~degree~C),
                                   title.position = "top", 
                                   title.hjust = 0.5,
                                   barwidth = unit(4, "in"), 
                                   frame.colour = "black", 
                                   ticks.colour = "black")) +
    map_theme +
    coord_sf(crs = world_moll) +
    labs(title =  "Global Shifts in Annual SST from Climatology Transition")


# Just plot them
stacked_clims <- (plot_82 + theme(legend.position = "none")) / plot_91
stacked_clims

plot_shift

Region Extents

The following figures look specifically at two regions. 1.) The Gulf of Maine, and 2.) the Northeastern United States’ Continental Shelf. Both of these areas have been mapped below for clarity on what data is included for any regional metric.

# load shapes/timeseries paths
gmri_areas    <- get_timeseries_paths(region_group = "gmri sst focal areas")
nelme_regions <- get_timeseries_paths(region_group = "nelme regions")

# get region shapes
gom          <- read_sf(gmri_areas$apershing_gulf_of_maine$shape_path)
ne_shelf     <- read_sf(nelme_regions$NELME$shape_path)



# lat/lon limits
map_lims <- list(x = c(-78, -63), y = c(35, 47))


# Gulf of Maine
gom_map <- base_map +
  geom_sf(data = gom, 
          alpha = 0.4,
          color = gmri_cols("gmri blue"),
          fill = gmri_cols("gmri blue")) +
  scale_x_continuous(breaks = seq(-180, 180, by = 5)) +
  coord_sf(xlim = map_lims$x, ylim = map_lims$y) +
  map_theme +
  labs(title = "Gulf of Maine")

#NE Shelf
shelf_map <- base_map +
  geom_sf(data = ne_shelf, alpha = 0.4,
          color = gmri_cols("teal"),
          fill = gmri_cols("teal")) +
  scale_x_continuous(breaks = seq(-180, 180, by = 5)) +
  coord_sf(xlim = map_lims$x, ylim = map_lims$y) +
  map_theme +
  labs(title = "Northeast Shelf")


# Mapping Them together
gom_map | shelf_map

Average Regional Temperatures

# masking function
mask_sf <- function(in_ras, sf_mask){
  r1 <- crop(in_ras, sf_mask)
  r2 <- mask(r1, sf_mask)
  return(r2)}

# Get Average Temp and Change in Climate
shift_table <- imap_dfr(grid_list, function(in_ras, climatology){
  # Get means
  gom_mean <- mask_sf(in_ras, gom) %>% cellStats(mean, na.rm = T)
  ne_mean  <- mask_sf(in_ras, ne_shelf) %>% cellStats(mean, na.rm = T)
  
  # build table
  grid_shift <- data.frame(
    Climatology = rep(climatology, 2),
    region = c("Gulf of Maine", "Northeast Shelf"),
    avg_temp = c(gom_mean, ne_mean))
  
  return(grid_shift) }) %>% 
    mutate(Climatology = factor(Climatology, levels = c("Climate Shift", "1991-2020", "1982-2011")))


# Plot the shifts
ggplot(shift_table, aes(x = avg_temp, y = Climatology)) +
  geom_segment(aes(xend = 0, yend = Climatology), color = gmri_cols("gmri blue")) +
  geom_point(size = 2, color = gmri_cols("orange")) +
  facet_wrap(~region) +
  geom_vline(xintercept = 0, color = "gray30", linetype = 3) +
  labs(y = "", x = expression("Average Temperature"~~degree~C)) 

Changes to Seasonality

Loading Climatology and Heatwave Info

# Load regional timeseries
gom_ts      <- read_csv(gmri_areas$apershing_gulf_of_maine$timeseries_path, col_types = cols()) %>% 
  mutate(time = as.Date(time))
ne_shelf_ts <- read_csv(nelme_regions$NELME$timeseries_path, col_types = cols()) %>%
  mutate(time = as.Date(time))


# Side by side - Heatwave Detection
# Run climatology for both periods, merge and compare
hw_comparison <- function(region_timeseries){
  
  # Run heatwave detection using new climate period
  heatwaves_82 <- pull_heatwave_events(region_timeseries, 
                                       threshold = 90, 
                                       clim_ref_period = c("1982-01-01", "2011-12-31")) %>% 
    select(time, sst, seas_82 = seas, anom_82 = sst_anom, hw_thresh_82 = mhw_thresh, cs_thresh_82 = mcs_thresh,
           hwe_82 = mhw_event, cse_82 = mcs_event)
           
  
  
  # Run heatwave detection using new climate period
  heatwaves_91 <- pull_heatwave_events(region_timeseries, 
                                       threshold = 90, 
                                       clim_ref_period = c("1991-01-01", "2020-12-31")) %>% 
    select(time, sst, seas_91 = seas, anom_91 = sst_anom, hw_thresh_91 = mhw_thresh, cs_thresh_91 = mcs_thresh,
           hwe_91 = mhw_event, cse_91 = mcs_event)
  
  
  # Subtract old from the new - set up flat date for plots
  base_date <- as.Date("2000-01-01")
  heatwaves_out <- left_join(heatwaves_82, heatwaves_91, by = c("time", "sst")) %>% 
    mutate(anom_shift = anom_91 - anom_82,
           clim_shift  = seas_91 - seas_82,
           upper_shift = hw_thresh_91 - hw_thresh_82,
           lower_shift = cs_thresh_91 - cs_thresh_82,
           hwe_change = case_when( 
             hwe_82 == T & hwe_91 == F ~ "No Longer HW",
             hwe_82 == F & hwe_91 == T ~ "New HW",
             hwe_82 == T & hwe_91 == T ~ "Still HW",
             hwe_82 == F & hwe_91 == F ~ "Still No HW",
             T ~ "missing case"),
           shift_direction = ifelse(clim_shift >= 0, "Temperature Increase", "Temperature Decrease"),
           year = year(time),
           yday = yday(time),
           flat_date = as.Date(yday-1, origin = base_date))
  
  
  return(heatwaves_out)
}


# Run comparisons
gom_comp <- hw_comparison(region_timeseries = gom_ts) 
ne_comp  <- hw_comparison(region_timeseries = ne_shelf_ts)

Plotting Seasonality

# Plotting the seasonal variations
seasonality_check <- function(comp_data){
  
  # Pull climatology
  seasonal_pattern <- comp_data %>% 
    filter(between(time, as.Date("2019-01-01"), as.Date("2019-12-31"))) %>% 
    distinct(flat_date, .keep_all = T)
  
  # Plot 82 clim
  p_82 <-ggplot(seasonal_pattern, aes(x = time)) + 
    geom_line(aes(y = seas_82), color = gmri_cols("gmri blue")) +
    geom_line(aes(y = hw_thresh_82), linetype = 2, color = "gray60") +
    geom_line(aes(y = cs_thresh_82), linetype = 2, color = "gray60") +
    labs(x = "", y = expression("SST"~~degree~C),
         subtitle = "1982-2011 Climatology") +
    scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0))
  
  # Plot 91 Clim
  p_91 <- ggplot(seasonal_pattern, aes(x = time)) + 
    geom_line(aes(y = seas_91), color = gmri_cols("teal")) +
    geom_line(aes(y = hw_thresh_91), linetype = 2, color = "gray60") +
    geom_line(aes(y = cs_thresh_91), linetype = 2, color = "gray60") +
    labs(x = "", y = expression("SST"~~degree~C),
         subtitle = "1991-2020 Climatology") +
    scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0))
  
  
  # Plot the difference
  p_diff <- ggplot(seasonal_pattern, aes(x = time, y = clim_shift)) + 
    geom_segment(aes(yend = 0, xend = time, color = shift_direction)) +
    scale_color_gmri(reverse = T) +
    labs(x = "", 
         y = expression("Temperature Shift"~~degree~C),
         color = "Shift in Climate Norm",
         subtitle = "Difference in Climatologies") +
    scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
    guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) +
    theme(legend.position = "bottom")
  
  return(
    list(
      "clim82" = p_82,
      "clim91" = p_91,
      "clim_shift" = p_diff)
  )
}

Gulf of Maine

# Return the plots
gom_seasonality <- seasonality_check(gom_comp)

# Assemble stack
(gom_seasonality[[1]] + labs(title = "Gulf of Maine")) / gom_seasonality[[2]] / gom_seasonality[[3]]

Northeast Shelf

# Return the plots
ne_seasonality <- seasonality_check(ne_comp)

# Assemble stack
(ne_seasonality[[1]] + labs(title = "Northeast Shelf")) / ne_seasonality[[2]] / ne_seasonality[[3]]

Heatwave Events

# Plotting how anomalies are now suppressed
heatwave_checks <- function(comp_data){
  
  # Set palette limits to center it on 0 with scale_fill_distiller
  limit <- max(abs(comp_data$anom_shift)) * c(-1,1)
  
  # Base date for rectangle fill for NA
  base_date <- as.Date("2000-01-01")
  
  # Assemble heatmap plot
  heatwave_heatmap <- ggplot(comp_data, aes(x = flat_date, y = year)) +
    
    # background box fill for missing dates
    geom_rect(xmin = base_date, 
              xmax = base_date + 365, 
              ymin = min(comp_data$year) - .5,
              ymax = max(comp_data$year) + .5, 
              fill = "gray75", 
              color = "transparent") +
    
    # Tile for sst colors
    geom_tile(aes(fill = anom_shift)) +
    
    # points for heatwave events
    geom_point(data = filter(comp_data, (hwe_82 == TRUE | hwe_91 == TRUE)),
               aes(x = flat_date, 
                   y = year,
                   color = hwe_change), 
               size = .25)  +
    
    scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
    scale_y_continuous(limits = c(1980.5, 2021.5), expand = c(0,0)) +
   
    
    # Colors
    scale_fill_gradient2(low =  gmri_cols("gmri blue"), high = "lightblue") + 
    scale_color_manual(values = c("black", "red")) + 
    
    #5 inches is default rmarkdown height for barheight
    guides("fill" = guide_colorbar(title = "Change in Anomaly Magnitude", 
                                   title.position = "top", 
                                   title.hjust = 0.5,
                                   barwidth = unit(3, "inches"), 
                                   frame.colour = "black", 
                                   ticks.colour = "black"),
           color = guide_legend(title = "Impact on Heatwave Record", 
                                title.position = "top", 
                                title.hjust = 0.5, 
                                override.aes = list(size = 3))) +  
     labs(x = "", 
          y = "") +
    theme_classic() +
    theme(legend.position = "bottom")


  # Assemble pieces
  return(heatwave_heatmap)
  
}

Gulf of Maine

heatwave_checks(gom_comp) + 
  labs(title = "Impact of New Climatology - Gulf of Maine")

Northeast Shelf

heatwave_checks(ne_comp) + 
  labs(title = "Impact of New Climatology - Northeast Shelf")

Net Drop of Heatwave Events